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Abstract 

We report on a novel stochastic analysis of seismic time series for the Earth's vertical velocity, by 
using methods originally developed for complex hierarchical systems, and in particular for turbulent 
flows. Analysis of the fluctuations of the detrended increments of the series reveals a pronounced 
change of the shapes of the probability density functions (PDF) of the series' increments. Before 
and close to an earthquake the shape of the PDF and the long-range correlation in the increments 
both manifest significant changes. For a moderate or large-size earthquake the typical time at which 
the PDF undergoes the transition from a Gaussian to a non-Gaussian is about 5-10 hours. Thus, 
the transition represents a new precursor for detecting such earthquakes. 

PACS number(s): 05.45.Tp 64.60.Ht 89.75.Da 91.30.Px 
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A grand challenge in geophysics, and the science of analyzing seismic activities, is devel- 
oping methods for predicting when earthquakes may occur. Any accurate information about 
an impending large earthquake may reduce the damage to the infrastructures and the number 
of casualties. Although there are many known precursors, experience over the last several 
decades indicates that reliable and quantitative methods for analyzing seismic data are still 
lacking [1]. 

A large number of stations around the world currently perform precise measurements 
of seismic time series. Several concepts and ideas have been advanced over the past few 
decades in order to explain important aspects of such seismic time series and what they imply 
for earthquakes, ranging from the Gutenberg-Richter law [2] for the number of earthquakes 
with a magnitude greater than a given value M, to the Omori law [3] for the distribution of 
the aftershocks, the concept of self-organized criticality [4], and the percolation model of the 
epicenters' spatial distribution [5]. At the same time, there has also been much interest in 
investigating the precursors to, and the predictability of, extreme increments in time scries 
[6] associated with disparate phenomena, ranging from earthquakes [7,8], to epileptic seizures 
[9], and stock market crashes [10-12]. 

In this Letter we provide compelling evidence for the existence of a novel transition in the 
probability density function (PDF) of the detrended increments of the stochastic fluctuations 
of the Earth's vertical velocity Vz, collected by broad-band stations (resolution 100 Hz). As 
an important new result, we demonstrate that there is a strong transition from a Gaussian 
to a non-Gaussian behavior of the increments' PDF as an earthquake is approached. We 
characterize the non-Gaussian nature of the PDF of the fluctuations in the increments of Vz, 
and the time- dependence of the PDF of the background fluctuations far from earthquakes. 

The results presented in this Letter are based on the detailed analysis of the data ob- 
tained from Spain's and California's broad-band networks for three earthquakes of large and 
intermediate magnitudes: the May 21, 2003, M = 7.1 event in Oran-Argel, detected in Ibiza 
(Balearic Islands); the 2004, M = 6.1 event in Alhucemas, and the M = 5.4 earthquake that 
occurred in California on April 30, 2008. Due to our recent discovery of localization of elastic 
waves in rock, both experimentally [13] and theoretically [14], we choose the distance of the 
detectors from the epicenters to be three times less than 300 km and one time 400 km. We 
have also analyzed many other earthquakes around the world, the results of which will be 
described briefly. 

The data are first detrended over different time scales in order to remove the possible 
trends in the time series x{t) = Vz{t). To do so, the time series is divided into semi-overlapping 
subintervals [1 + s{k — 1), s{k + 1)] of length 2s and labeled by the index k > 1. Next, we 
fit x{t) to a third-order polynomial [15-17] to detrend the original series in the corresponding 
time window. The detrended increments on scale s are defined by Zs{t) = x*{t + s) — x*{t), 
where t E [1 + s{k — 1), sk], with x*{t) being the detrended series, i.e., the deviation of x{t) 
from its fitted value. 

We then develop a new approach, originally proposed for fully-developed turbulence 
[18-21], in order to describe the cascading process that determines how the fluctuations in the 
series evolve, as one passes from the coarse to fine scales. For a fixed t, the fiuctuations at 
scales s and As are related through the cascading rule, 

Zxs{t) = WxZsit), Vs, A > , (1) 
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where In(WA) is a random variable. Iterating Eq. (1) forces implicitly the random variable Wx 
to follow a log infinitely-divisible law [22] . One of the simplest candidates for such processes 
is represented [17] by, Zs{t) — Cs(^) exp[a;s(f)], where and u!s{t) are independent Gaussian 
variables with zero mean and variances and a^. The PDF of Zs{t) has fat tails that depend 
on the variance of cUs, and is expressed by [21], 

Ps{Z,) ^ j ^Gs{\na)d\na , (2) 

where Fg and Gg are both Gaussian with zero mean and variances and A^, respectively, 
e.g., ^^(Ina-) = l/(-\/27rAs) exp(— In^ (7/2A^). In this case, Ps{Zs) is expressed by Eq. (2) and 
converges to a Gaussian distribution as A^ — > 0. Although Eq. (2) is equivalent to that for a 
log-normal cascade model, originally introdTiccd to study fully-developed turbulence [23,24], 
it also describes approximately the non-Gaussian PDFs observed in a broad range of other 
phenomena and systems, such as foreign exchange markets [17,24,25] and heartbeat interval 
fluctuations [17,18] (see also [26-28]). 

To carry out a quantitative analysis of the seismic times series, we focus on two aspects, 
namely, deviations of the PDF of the detrended increments from a Gaussian distribution, and 
the dependence of correlations in the increments on the scale parameter s. We begin with 
the time series for the largest (Af = 7.1) earthquake for two distinct time intervals: (i) data 
set (I) representing the background fluctuations far from the time of the earthquake, and (ii) 
data set (II) close (less than 5 hours) to the earthquake. 

To flt the increments' PDF to Eq. (2), wc estimate the variance \1{s), using the least- 
squares method, with the error bars estimated by the goodness of the fit method. Deviation 
of A^(s) from zero is a possible indicator of non-Gaussian statistics. As shown in Fig. 1, we 
find an accurate parametrization of the PDFs by \1{s) for both data sets. Moreover, the PDF 
of Zg for the data set (I) becomes essentially Gaussian as s increases to 800 ms, whereas it 
deviates from the Gaussian distribution for the data set (II). The time scale s = 800 ms for A^ 
within a moving window was estimated by the plot of A^ vs s for the data set (I) (background 
fiuctuations), and selecting s such that A^ (see also below). 

The scale-dependence of the parameter A^ of the PDF is shown in Fig. 2. For the data 
set (I) of the M — 7.1 earthquake, shown in Fig. 2(a), and times 200 ms < s < 500 ms, 
a logarithmic decay, A^ oc logs, is obtained. For the data set (II), the logarithmic regime 
extends to 300 ms < s < 2000 ms. Figure 2(b) presents similar behavior for the M — 6.1 
earthquake. We note that, for the data set (II) of the M — 6.1 earthquake, there is a crossover 
time at which A^ changes from a ~ log(s) behavior to having a finite value, ~ 0.3. 

The importance of the results shown in Fig. 2 is that, they indicate that the increments' 
PDFs for s > 2000 ms and s > 1500 ms are almost Gaussian (A^ — > 0) for the M = 7.1 and 
M — 6.1 earthquakes [for the data set (II)], respectively. Transforming the time scales to 
length scales via the velocity of the elastic waves in Earth, ~ 5000 m/sec, the corresponding 
length scales are about 10 km and 7.5 km, for the same earthquakes, respectively, implying 
that larger earthquakes have larger characteristic length scales, and that for the M — 6.1 
event, there is a smaller active part in the fault. 

We note that as one moves down the cascade process from the large to small scales, one 
expects the statistics to increasingly deviate from Gaussianity (see above and Fig. 2), in order 
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to derive Eq. (2). Note that, a non-Gaussian PDF with fat tails on small scales indicates an 
increased probability of the occurrence of short-time extreme seismic fluctuations. 

Prom the point of view of the increments' PDF, the non-Gaussian noise with uncorre- 
lated oUg in the process, Zs{t) — ^^(i) exp[a;s(t)], and a multifractal formulation are indistin- 
guishable, because their one-point statistics at any given scale may be identical. Thus, to 
understand the origin of the non-Gaussian fluctuations, we explore the correlation properties 
of uig. To do so, we use an alternative method for studying the correlation functions of the 
local "energy" fluctuations [20] . We define the magnitude of local variance over a scale s by, 
^si^)s^ ^'k=-ns/2 ^s{t + kAty, and, a)s(i) = | log cr^ (i), respectively. Here, At is the sampling 
interval and, Ug = s/ At. The magnitude of the correlation function of cUg is then defined by 

C^'\t) = {[us{t) - {ujs)][u,s{t + r) - {ujs)]) , (3) 

where (•) indicates a statistical average. Figure 3 shows the results for the two data sets. The 
correlation function decays sharply for the data set (I) - far from the earthquakes - whereas 
it is of long-range type for the set (11) - close to the earthquakes. 

We emphasize that for the data set (11), the PDF deviates from being Gaussian even 
for s > 800 ms. Although one might argue that the deviations might be due to an underlying 
Levy statistics, this possibility is ruled out due to the deduced hierarchical structures that 
imply that the increments for different scales are not independent; see Fig 3. 

We now show that the analysis may be used as a new precursor for detecting an impend- 
ing earthquake. A window containing one hour of data is selected and moved with At = 15 
minutes to determine the temporal dependence of A^. Guided by Fig. 2, the local temporal 
variations of for s = 800 ms are investigated. According to Fig. 2, for s ^ 800 ms, the 
difference between the values of A^ is large enough for the background data and the data set 
near the earthquakes. Hence, such a time scale may be used as the characteristic time for the 
dynamics of the non-Gaussian indicator A^. Figures 4(a) and 4(b) display a well- pronounced, 
systematic increase in A^ as the earthquakes are approached. Taking into account the esti- 
mated error of A^ for the background fiuctuations in Figs. 4, wc sec that about 7 and 5 hours 
before the earthquakes values of A^ are larger, by more than two standard deviations, than 
those for the background. 

We note that in defining A^, one supposes that the PDF of the increments is log- normal. 
This may induce errors due to the fact that, one must fit the PDF with a predefined functional 
form. An unbiased quantity that measures the intermittency and deviation from Gaussian is 
the fiatness, which provides an unbiased estimator of deviations from Gaussianity, without 
having to adopt a priori any functional form for the PDF. In Figs. 4(c) and 4(d) we present 
the flatness of the time series in the same windows (see above) for time scale s ~ 800 ms. 
They show that, consistent with A^, the flatness also yields a clear alert for an impending 
earthquake. 

Due to the localization of elastic waves in Earth, stations that are far from an earthquake 
epicenter cannot provide any clue to the transition in the shape of the increments' PDF. We 
checked this point for several earthquakes. Shown in Fig. 5 are the results for the M — 5.4 
California earthquake, occurred at (40.837 N, 123.499 W). The distance from the epicenter of 
the station that does provide an alert of about 3 hours for the earthquake is about ~ 128 km, 
while the second station that does not yield any alert is about ~ 400 km from the epicenter. 
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We also analyzed several other earthquakes of various magnitudes. We found that for 
earthquakes with M < 5 the increase in is not large, even for the data that are collected 
in stations about 100 km from the epicenters. Moreover, we also analyzed seismic data for 
large earthquakes in Pakistan and Iran, and found that they exhibit the same types of trends 
and results, as those presented above, for the time variation of close to the earthquakes. 
For example, for the M = 7.6 earthquake that occurred on August 10, 2005, in Pakistan, the 
transition in the value of A^ occurred about 10 hours before the earthquake, while for the 
M — 6.3 earthquake that occurred in northern Iran on May 28, 2004, the transition happened 
about 4 hours before the earthquake. 

In summary, the temporal dependence of the fat tails of the PDF of the increments of 
the vertical velocity Vz{t) of Earth exhibits a gradual, systematic increase in the probability 
of the appearance of large values on approaching a large or moderate earthquake, which is 
interpreted as an alert for the earthquake. To estimate the alert time one must, (i) utilize the 
time scries Vz{t), collected at broad-band stations near the epicenters. Due to localization of 
elastic waves in rock [15], data from far away stations cannot provide any clue to the transition 
in the shape of the increments' PDF. The station's distance (300 km) is not universal and 
depends on the geology, but it is of the correct order of magnitude, (ii) The time scale s for 
moving A^ (here, 800 ms) is estimated by its plot vs s for data set (I), and a choice of s such 
that A^ — > 0. On this scale the difference between A^ for the data sets (I) and (II) will be 
large enough to obtain a meaningful alert for the earthquake, (iii) One must also estimate A^ 
or the flatness in some windows (here, 1 hour windows) and move it over the time series, in 
order to observe its variations with the time. 
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FIG. 1: Continuous deformation of the increments' PDFs in log-linear scale, for the M = 7.1 
earthquake across scales for, from top to bottom, s = 200, 400, 600, and 800 ms, and (a) far from, 
and (b) close to, the earthquake. Solid curves are the PDFs based on Eq. (2), while dashed curves 
are the Gaussian PDF. 
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FIG. 2: Scale-dependence of Ag vs logs, (a) The M = 7.1 event, both far from [data set (I)] and 
close to [data set (II)] the earthquake. For the data set (I) and s > 700 ms, — > 0, implying that 
the increments' PDF is Gaussian, whereas for the data set (II) A^ deviates strongly from zero for 
700 ms < s < 1500 ms. (b) Same as in (a), but for the M = 6.1 earthquake. When A^ — > the error 
bars are very small, about the same size as the symbols. 
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FIG. 3: The correlation function C^^\t) for the data set (I) far from, and the set (II) close to, the 
M = 7.1 earthquake. 
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FIG. 4: The local temporal dependance of and the flatness for s = 800 ms, over a one-hour period, 
for the M = 7.1 and M = 6.1 events, indicating a gradual, systematic increase on approaching the 
earthquakes. 
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FIG. 5: (a) The data for the M = 5.4 earthquake in Cahfornia. (b) and (c) Show the local temporal 
dependance of for s = 500 ms, collected at stations with a distance d from the epicenter. The 
station at d = 128 km provides the alert, whereas the second station does not yield any alert. 
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